The first thing to note is that, because the data was from a variety of data sources from the UK government then they were in different geographies as well. While those that were set in 2018 or 2016 LADs were not a problem to convert to NUTS3 geography as the lookups were readily available, the indices of multiple deprivation and unemplyment data were set in 2019 LAD geography which as of yet did not have a readily available lookup. Therefore, to map the unemployment and the indices of multiple derpivation onto the NUTS3 region I had to do a point in polygon analysis to find the centroids of the LAD19 boundaries and see what NUTS3 region they were located within. Once this was done they could then be aggregated up to the NUTS3 geography along with the indicators. Below is the code that was used to find these geographies.

library(sf)
## Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3
library(tmap)
library(leafpop)
library(leaflet)
library(tmaptools)
library(tidyverse)
## -- Attaching packages ----------------- tidyverse 1.2.1 --
## v ggplot2 3.2.1     v purrr   0.3.2
## v tibble  2.1.3     v dplyr   0.8.3
## v tidyr   1.0.0     v stringr 1.4.0
## v readr   1.3.1     v forcats 0.4.0
## -- Conflicts -------------------- tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(plyr)
## -------------------------------------------------------------------------
## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)
## -------------------------------------------------------------------------
## 
## Attaching package: 'plyr'
## The following objects are masked from 'package:dplyr':
## 
##     arrange, count, desc, failwith, id, mutate, rename, summarise,
##     summarize
## The following object is masked from 'package:purrr':
## 
##     compact
library(classInt)
library(RColorBrewer)
library(geojsonio)
## Registered S3 method overwritten by 'geojsonio':
##   method         from 
##   print.location dplyr
## 
## Attaching package: 'geojsonio'
## The following object is masked from 'package:base':
## 
##     pretty
library(knitr)

#importing the NUTS3 geojson
NUTS3 <- geojson_read("https://opendata.arcgis.com/datasets/473aefdcee19418da7e5dbfdeacf7b90_2.geojson", what = "sp")

#switching the geojson to NUTS
NUTS3_SF <- st_as_sf(NUTS3)

#The 2019 LAD do not have projections so I must attempt to draw LAD19s onto NUTS3 projectsion
#To do this do a point in polygon analysis
#get the LAD19 gejson
LAD19 <- geojson_read("https://opendata.arcgis.com/datasets/cec4f9cf783a47bab9295b2e513dd342_0.geojson", what = "sp")
#https://opendata.arcgis.com/datasets/2b95585accc4437b97d766f31c5568cb_0.geojson

#shift this to SF
LAD19_SF <- st_as_sf(LAD19)

#make sure both the NUTS3 map and LAD19 map have the same projection
LAD19_SF <- st_transform(LAD19_SF, 27700)
NUTS3_SF <- st_transform(NUTS3_SF, 27700)

#convert the LAD19 boundaries into points
LAD19_centroids <- st_centroid(LAD19_SF)
## Warning in st_centroid.sf(LAD19_SF): st_centroid assumes attributes are
## constant over geometries of x
#set the CRS for both of these
st_crs(LAD19_centroids) <- 27700
st_crs(NUTS3_SF) <- 27700

#join them on the line at which they intersect
Joined <- st_join(LAD19_centroids, NUTS3_SF, join = st_intersects)

#create a df of the joined sf df so that I can use VLOOKUP in excel to map LAD19s onto 
df <- Joined[,c("LAD19NM", "LAD19CD", "LAD19NMW", "nuts318cd")]


#remove the geometry from the data set
df_df <- df %>% st_set_geometry(NULL)
#check that it has been removed
class(df_df)
## [1] "data.frame"
#write the dataset to a csv to be able to use VLOOKUP in excel to map LAD19s onto NUTS3 regions
write.csv(df_df, 'LAD19.csv')

#the mapping of boundaries of one onto another was done outside of R and in excel
#this was due to the difficulty of doing so in R
#They were mapped using the Countif and sumif functions to average out if moving up into a higher mapping level

Once this was done and all the measures were harmonised in terms of geography they were then inputted into one spreadsheet. This spreadsheet was then imported into R using the read.csv function and setting and missing values as na. These were then merged to create the NUTS3Map

#reading in the data that was created as a result
NUTS3_data <- read.csv("Data/NUTS3(1).csv", na=c("n/a", "na"))

#merging the sf object with the data created to be able to analyse this later
UKNUTS3Map <- merge(NUTS3_SF,
                    NUTS3_data,
                    by.x = "nuts318cd",
                    by.y = "NUTS.code",
                    no.dups=TRUE)

The next thing to do was to get the breaks for the plots. This was done by finding the means, standard deviations, minimum and maximum for the relevant measures as seen below.

#setting the breaks for the GVA data
breaksGVA <- c(0,55,70,85,100,115,130,145,1500)

breaksGVA
## [1]    0   55   70   85  100  115  130  145 1500
#finding the mean and standard deviation for the relevant columns to be able to set the breaks
UEM <- mean(UKNUTS3Map$U.E.rate...Jul.2018.Jun.2019, na.rm=TRUE)
UESD <- sd(UKNUTS3Map$U.E.rate...Jul.2018.Jun.2019, na.rm=TRUE)
UEMin <- min(UKNUTS3Map$U.E.rate...Jul.2018.Jun.2019, na.rm=TRUE)
UEMax <- max(UKNUTS3Map$U.E.rate...Jul.2018.Jun.201, na.rm=TRUE)

#the breaks were then set using the standard deviations and means
breaksUE <- c(UEMin, UEM-2*UESD, UEM-1.25*UESD, UEM-0.5*UESD, UEM, UEM+0.5*UESD, UEM+1.25*UESD, UEM+2*UESD, UEMax)

breaksUE
## [1] 0.800000 1.740037 2.636242 3.532448 4.129918 4.727389 5.623594 6.519800
## [9] 8.100000
#This was then repeated for the other measures e.g. IMD
IMDM <- mean(UKNUTS3Map$IMD.2019...Average.score, na.rm=TRUE)
IMDSD <- sd(UKNUTS3Map$IMD.2019...Average.score, na.rm=TRUE)
IMDMax <- max(UKNUTS3Map$IMD.2019...Average.score, na.rm=TRUE)
IMDMin <- min(UKNUTS3Map$IMD.2019...Average.score, na.rm=TRUE)

breaksIMD <- c(IMDMin, IMDM-1.5*IMDSD, IMDM-0.75*IMDSD, IMDM, IMDM+0.75*IMDSD, IMDM+1.5*IMDSD, IMDMax)

breaksIMD
## [1]  9.58950 10.86407 16.53617 22.20828 27.88038 33.55249 45.03900
#male life expectancy
LEMM <- mean(UKNUTS3Map$LE.males.2015.2017,na.rm=TRUE)
LEMSD <- sd(UKNUTS3Map$LE.males.2015.2017, na.rm=TRUE)
LEMMin <- min(UKNUTS3Map$LE.males.2015.2017, na.rm=TRUE)
LEMMax <- max(UKNUTS3Map$LE.males.2015.2017, na.rm=TRUE)

breaksLEM <- c(LEMMin, LEMM-1.5*LEMSD, LEMM-0.75*LEMSD, LEMM, LEMM+0.75*LEMSD, LEMM+1.5*LEMSD, LEMMax)

breaksLEM
## [1] 74.20000 77.23186 78.29379 79.35573 80.41766 81.47959 82.30000
#female life expectancy
LEFM <- mean(UKNUTS3Map$LE.females.2015.2017, na.rm = TRUE)
LEFSD <- sd(UKNUTS3Map$LE.females.2015.2017, na.rm = TRUE)
LEFMin <- min(UKNUTS3Map$LE.females.2015.2017, na.rm=TRUE)
LEFMax <- max(UKNUTS3Map$LE.females.2015.2017, na.rm=TRUE)

breaksLEFM <- c(LEFMin, LEFM-1.5*LEFSD, LEFM-0.75*LEFSD, LEFM, LEFM+0.75*LEFSD, LEFM+1.5*LEFSD, LEFMax)

breaksLEFM
## [1] 79.51983 81.03161 82.03371 83.03582 84.03792 85.04002 86.54595
#GCSE scores
GCSEM <- mean(UKNUTS3Map$GCSE.2018.A..C.., na.rm=TRUE)
GCSESD <- sd(UKNUTS3Map$GCSE.2018.A..C.., na.rm=TRUE)
GCSEMin <- min(UKNUTS3Map$GCSE.2018.A..C.., na.rm=TRUE)
GCSEMax <- max(UKNUTS3Map$GCSE.2018.A..C.., na.rm=TRUE)

breaksGCSE <- c(GCSEMin, GCSEM-2*GCSESD, GCSEM-GCSESD, GCSEM, GCSEM+GCSESD, GCSEM+2*GCSESD, GCSEMax)

breaksGCSE
## [1] 45.60000 52.84907 58.18581 63.52256 68.85930 74.19605 76.70000
#percentage voting leave
LeaveM <- mean(UKNUTS3Map$Leave, na.rm=TRUE)
LeaveSD <- sd(UKNUTS3Map$Leave, na.rm=TRUE)
LeaveMin <- min(UKNUTS3Map$Leave, na.rm=TRUE)
LeaveMax <- max(UKNUTS3Map$Leave, na.rm=TRUE)

breaksLeave <- c(LeaveMin, LeaveM-1.5*LeaveSD, LeaveM-0.75*LeaveSD, LeaveM, LeaveM+0.75*LeaveSD, LeaveM+1.5*LeaveSD, LeaveMax)

breaksLeave
## [1] 21.40000 38.20013 45.88428 53.56842 61.25257 68.93671 72.30000
#based on these brakes we could then get the colour palettes that would be associated with the relevant plots
RdBu7 <- get_brewer_pal("RdBu", n=7)

RdBu5 <- get_brewer_pal("RdBu", n=5)

RdBu8 <- get_brewer_pal("RdBu", n=8)

ReverseRdBu <- get_brewer_pal("-RdBu", n =6)

ReversedRdBU <- get_brewer_pal("-RdBu", n=8)

Based on these breaks we could then draw the following plots. The first to be drawn is that of GVA which was used to inform the third and final boundary measure as seen below.

tm1 <- tm_shape(UKNUTS3Map)+
  #setting the fill as GVA
  tm_fill("GVA.in..2017",
          #setting the alpha to 0
          alpha = 1,
          #setting the title of the legend
          title= "GVA (% of UK average)",
          #setting the breaks
          breaks=breaksGVA, 
          #selecting the pallette
          palette = RdBu8,
          #setting the legend to show
          legend.show = TRUE)+
  #setting borders between the regions
  tm_borders(
             #setting the colour as grey 
             col = "grey",
             #setting the alpha as 0.8 
             alpha= 0.8)+
  #adding a legend
  tm_legend(#setting the title size
            title.size = 1.2,
            show=TRUE, 
            #setting the position
            position = c(0.05,0.4))+
  tm_layout(#setting the title for the overall plot
            title = "(a)",
            #setting the size of the title
            title.size = 2,
            #setting the title's position
            title.position = c("left", "top"),
            #setting a frame around the plot
            frame=TRUE,
            #setting the heigh of the legend
            legend.height = 4,
            #setting its width
            legend.width= 4,
            #setting the text size for the legend as 1
            legend.text.size = 1)+
  #adding a compass to the plot
  tm_compass(#choosing the type of the compass
             type = "rose",
             #positioning the compassin the top right color
             position = c(0.75, 0.8,
             #setting the color
             color.light="grey90"),
             #the usual size is 6 but we want it smaller
             size = 4)+
  #adding a scale bar to the plot
  tm_scale_bar(#setting the breaks as 0,50 and 100km
               breaks = c(0,50,100),
               text.size=1)

tm1

We are then concerned with drawing the different boundary lines. Following inspection of Rowthron’s (2010) paper and Dorling’s (2010) paper two boundaries can be drawn. The first uses the boundaries of the NUTS1 regions in the UK, the second uses Parliamentary constituencies. While the first can be mapped exactly, the second must be mapped approximately based on how the outline fits with the NUTS3 boundaries. The third plot is then deduced from a visual inspection of the previous plot, along with ideas from the national media suggesting that the South-East is seperating from the rest of the UK.

tm2 <- tm_shape(UKNUTS3Map)+
  #setting a fill along the dividing line, 1 = North, 0 = South
  tm_fill("Dorling")+
  #showing that we don't want a legend
  tm_legend(show=FALSE)+
  #showing that we don't want a frame
  tm_layout(frame=FALSE)+
  #setting the credit as a)
  tm_credits("(b)",
             #setting the position and size of this credit
             position = c(0.3,0.80),
             size=1.5)

#Repeating this for Rowthorn and South-East definitions
tm3 <- tm_shape(UKNUTS3Map)+
  tm_fill("Rowthorn")+
  tm_legend(show=FALSE)+
  tm_layout(frame=FALSE)+
  tm_credits("(c)",
             position = c(0.3,0.80),
             size=1.5)

tm4 <- tm_shape(UKNUTS3Map)+
  tm_fill("SE")+
  tm_legend(show=FALSE)+
  tm_layout(frame=FALSE)+
  tm_credits("(d)", 
             position = c(0.3,0.80), 
             size=1.5)

#arranging these plots in a single column so they can be easily compared
tmap_arrange(tm2,tm3,tm4, ncol=1)

This can then be followed up by plots for the other six indicators

#creating a plot to show unemployment
tm5 <- tm_shape(UKNUTS3Map)+
  #creating a fill
  tm_fill(#setting where to get the data from
          "U.E.rate...Jul.2018.Jun.2019",
          #setting the alpha to 1
          alpha = 1,
          #setting the title of the legend
          title= "Unemployment (%)",
          #setting the breaks as what has been created before
          breaks= breaksUE, 
          #setting the pallete as chosen before
          palette = ReverseRdBu,
          legend.width = 4,
          #showing the legend
          legend.show = TRUE)+
  #creating borders 
  tm_borders(#setting the colour as grey
             col = "grey",
             #setting the alpha as 0.5 so it isn't as strong
             alpha= 0.5)+
  #adding the legend
  tm_legend(#setting the title size to 1
            title.size = 1,
            show=TRUE, 
            #setting the position of the title
            position = c(0.05,0.3)
  )+
  tm_layout(#setting the overall title
            title = "Unemployment Rate 2018-2019",
            #selecting the title position as the left and top
            title.position = c("left", "top"),
            #selecting that we only want 1 digit after the decimal place
            legend.format=list(fun = function(x) paste0(formatC(x, digits = 1, format = "f"))),
            #selectring that we want a frame for the overall plot
            frame=TRUE,
            legend.title.size = 1,
            legend.show = TRUE,
            legend.position = c(0.05, 0.35),
            legend.width=3)+
  #adding a compass
  tm_compass(#setting the type as a rose
             type = "rose",
             #setting the position in the top right corner
             position = c(0.75, 0.8,
                          color.light="grey90"),
             #setting the size as 3
             size = 3)+
  #adding a scale bar
  tm_scale_bar(#setting the breaks as 0,50 and 100
               breaks = c(0,50,100),
               #setting the text size as 1
               text.size=1)

tm5

tm10 <- tm_shape(UKNUTS3Map)+
  tm_fill("IMD.2019...Average.score",
          alpha = 1,
          title= "IMD score",
          breaks= breaksIMD, 
          palette = ReverseRdBu,
          legend.width = 4,
          legend.show = TRUE)+
  tm_borders(col = "grey",
             alpha= 0.5)+
  tm_legend(title.size = 1,
            show=TRUE, 
            position = c(0.05,0.35))+
  tm_layout(title = "IMD 2019 average score",
            title.position = c("left", "top"),
            legend.format=list(fun = function(x) paste0(formatC(x, digits = 1, format = "f"))),
            frame=TRUE)+
  tm_compass(type = "rose",
             position = c(0.75, 0.8,
                          color.light="grey90"),
             size = 3)+
  tm_scale_bar(breaks = c(0,50,100),
               text.size=1)

breaksIMD
## [1]  9.58950 10.86407 16.53617 22.20828 27.88038 33.55249 45.03900
tm10

tm6 <- tm_shape(UKNUTS3Map)+
  tm_fill("LE.females.2015.2017",
          alpha = 1,
          title= "Life expactancy (years)",
          breaks= breaksLEFM, 
          palette = RdBu7,
          legend.width = 4,
          legend.show = TRUE)+
  tm_borders(col = "grey",
             alpha= 0.5)+
  tm_legend(title.size = 1,
            show=TRUE, 
            position = c(0.05,0.4))+
  tm_layout(title = "Female life expectancy 2015-2017 ",
            title.position = c("left", "top"),
            legend.format=list(fun = function(x) paste0(formatC(x, digits = 1, format = "f"))),
            frame=TRUE)+
  tm_compass(type = "rose",
             position = c(0.75, 0.8,
                          color.light="grey90"),
             size = 3)+
  tm_scale_bar(breaks = c(0,50,100),
               text.size=1)

tm6

tm7 <- tm_shape(UKNUTS3Map)+
  tm_fill("LE.males.2015.2017",
          alpha = 1,
          title= "Life expectancy (years)",
          breaks= breaksLEM, 
          palette = RdBu7,
          legend.width = 4,
          legend.show = TRUE)+
  tm_borders(col = "grey",
             alpha= 0.5)+
  tm_borders(col = "grey",
             alpha= 0.5)+
  tm_legend(title.size = 1,
            show=TRUE, 
            position = c(0.05,0.35))+
  tm_layout(title = "Male Life expactancy 2015-2017",
            title.position = c("left", "top"),
            legend.format=list(fun = function(x) paste0(formatC(x, digits = 1, format = "f"))),
            frame=TRUE)+
  tm_compass(type = "rose",
             position = c(0.75, 0.8,
                          color.light="grey90"),
             size = 3)+
  tm_scale_bar(breaks = c(0,50,100),
               text.size=1)

tm7
## Warning: One tm layer group has duplicated layer types, which are omitted.
## To draw multiple layers of the same type, use multiple layer groups (i.e.
## specify tm_shape prior to each of them).

#hist(NUTS3_data$GCSE.2018.A..C.., main="Histogram of unemployment", xlab="Unemployment", xlim=c(45,78), breaks=30)


tm8 <- tm_shape(UKNUTS3Map)+
  tm_fill("GCSE.2011.A..C..",
          alpha = 1,
          title= "A*-C grades (%)",
          breaks= breaksGCSE, 
          palette = RdBu7,
          legend.width = 4,
          legend.show = TRUE)+
  tm_borders(col = "grey",
             alpha= 0.5)+
  tm_legend(title.size = 1,
            show=TRUE, 
            position = c(0.05,0.4))+
  tm_layout(title = "Percentage of A*-C GCSE grades 2018",
            title.position = c("left", "top"),
            legend.format=list(fun = function(x) paste0(formatC(x, digits = 1, format = "f"))),
            frame=TRUE)+
  tm_compass(type = "rose",
             position = c(0.75, 0.8,
                          color.light="grey90"),
             size = 3)+
  tm_scale_bar(breaks = c(0,50,100),
               text.size=1)

tm8
## Warning: Values have found that are less than the lowest break

#hist(NUTS3_data$Leave, main="Histogram of unemployment", xlab="Unemployment", xlim=c(20,80), breaks=100)

tm9 <- tm_shape(UKNUTS3Map)+
  tm_fill("Leave",
          alpha = 1,
          title= "Leave (%)",
          breaks= breaksLeave, 
          palette = ReverseRdBu,
          legend.width = 4,
          legend.digits = 3,
          legend.show = TRUE)+
  tm_borders(col = "grey",
             alpha= 0.5)+
  tm_legend(title.size = 1,
            show=TRUE, 
            position = c(0.05,0.4))+
  tm_layout(title = "Percentage voting Leave in 2016",
            title.position = c("left", "top"),
            legend.format=list(fun = function(x) paste0(formatC(x, digits = 1, format = "f"))),
            frame=TRUE)+
  tm_compass(type = "rose",
             position = c(0.75, 0.8,
                          color.light="grey90"),
             size = 3)+
  tm_scale_bar(breaks = c(0,50,100),
               text.size=1)

tm9

#these six plots are then arranged into a single plot such that there are three columns with 2 layers
t = tmap_arrange(tm5, tm6, tm7,tm8,tm10,tm9, ncol=3)
t
## Legend labels were too wide. The labels have been resized to 0.66, 0.66, 0.66, 0.66, 0.66, 0.66. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.
## Warning: One tm layer group has duplicated layer types, which are omitted.
## To draw multiple layers of the same type, use multiple layer groups (i.e.
## specify tm_shape prior to each of them).
## Some legend labels were too wide. These labels have been resized to 0.66, 0.66, 0.66, 0.66, 0.66, 0.66. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.
## Warning: Values have found that are less than the lowest break
## Legend labels were too wide. The labels have been resized to 0.66, 0.66, 0.66, 0.66, 0.66, 0.66. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.
## Some legend labels were too wide. These labels have been resized to 0.66, 0.66, 0.66, 0.66, 0.66. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.
## Legend labels were too wide. The labels have been resized to 0.66, 0.66, 0.66, 0.66, 0.66, 0.66. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.

#this ouput is then exported to produce a png

These static plots are included in the PDF output however these can also be mapped in an interactive plot that includes all seven indicators.

library(leafpop)
library(leaflet)

#given that leaflet uses world map data then the map must be ttrasnsformed to WGS84 projection
UKNUTS3MapWGS <- st_transform(UKNUTS3Map, 4326)

#the data can be used to create popup tables for each of the indicators
popGVA <- popupTable(#the data comes from the map
                     UKNUTS3MapWGS,
                     #the columns to be included are the NUTS3 code, the NUTS3 name and the valye
                     zcol=c("nuts318cd", "nuts318nm", "GVA.in..2017"))
#unemployment popup
popUE <- popupTable(UKNUTS3MapWGS, zcol=c("nuts318cd", "nuts318nm", "U.E.rate...Jul.2018.Jun.2019"))
#male life expectancy popup
popMLE <- popupTable(UKNUTS3MapWGS, zcol=c("nuts318cd", "nuts318nm", "LE.males.2015.2017"))
#female life expectancy popup
popFLE <- popupTable(UKNUTS3MapWGS, zcol=c("nuts318cd", "nuts318nm", "LE.females.2015.2017"))
#GCSE popup
popGCSE <- popupTable(UKNUTS3MapWGS, zcol=c("nuts318cd", "nuts318nm", "GCSE.2018.A..C.."))
#IMD popup
popIMD <- popupTable(UKNUTS3MapWGS, zcol=c("nuts318cd", "nuts318nm", "IMD.2019...Average.score"))
#Brexit popup
popBrexit <- popupTable(UKNUTS3MapWGS, zcol=c("nuts318cd", "nuts318nm", "Leave"))

#The data is then used to set the colour pallets to be used
#palette 1 for GVA
pal1 <- colorBin(palette="RdBu",
                 #where the data comes from
                 domain=UKNUTS3MapWGS$GVA.in..2017,
                 #what breaks to use
                 bins=breaksGVA)
#palette for unemployment
pal2 <- colorBin(palette=ReverseRdBu, domain=as.numeric(UKNUTS3MapWGS$U.E.rate...Jul.2018.Jun.2019), bins=breaksUE)
#palette for IMD
pal3 <- colorBin(palette=ReverseRdBu, domain=as.numeric(UKNUTS3MapWGS$IMD.2019...Average.score), bins=breaksIMD)
#palette for female life expectancy
pal4 <- colorBin(palette="RdBu", domain=as.numeric(UKNUTS3MapWGS$LE.females.2015.2017), bins=breaksLEFM)
#palette for male life expectancy
pal5 <- colorBin(palette="RdBu", domain=as.numeric(UKNUTS3MapWGS$LE.males.2015.2017), bins=breaksLEM)
#palette for GCSE scores
pal6 <- colorBin(palette = "RdBu", domain=as.numeric(UKNUTS3MapWGS$GCSE.2018.A..C..), bins = breaksGCSE)
#palette for the percentage voting leave
pal7 <- colorBin(palette = ReverseRdBu, domain=as.numeric(UKNUTS3MapWGS$Leave), bins = breaksLeave)

#these popup tables and palettes can then be called in the leaflet map

#creating the leaflet map
map <- leaflet(UKNUTS3MapWGS) %>%
  
  #creating basemap options
  addTiles(group = "OSM (default)") %>%
  #adding polygons
  addPolygons(
              #fillcolor comes from the palette defined before and the data is GVA
              fillColor = ~pal1(GVA.in..2017),
              #setting the base color
              color = "white",
              weight = 2,
              opacity = 1,
              dashArray = "3",
              #setting the popup table defined before
              popup = popGVA,
              #setting how opaque the fills are
              fillOpacity = 0.7,
              #setting the group to link the legend
              group = "GVA",
              #adding a highlight option
              highlight = highlightOptions(
                dashArray = "",
                #setting that it becomes more noticeable
                fillOpacity = 0.8,
                weight = 3,
                #setting the colour of the outline
                color = "Grey",
                #bringing the plot to the front
                bringToFront = TRUE)
              ) %>%
  
  #this is then replicated for the other indicators
  addPolygons(fillColor = ~pal2(as.numeric(U.E.rate...Jul.2018.Jun.2019)),
              color = "white",
              weight = 2,
              opacity = 1,
              dashArray = "3",
              popup = popUE,
              fillOpacity = 0.7,
              group = "UE",
              highlight = highlightOptions(
                dashArray = "",
                fillOpacity = 0.8,
                weight = 3,
                color = "Grey",
                bringToFront = TRUE)) %>%
  
  addPolygons(fillColor = ~pal3(IMD.2019...Average.score),
              color = "white",
              weight = 2,
              opacity = 1,
              dashArray = "3",
              popup = popIMD,
              fillOpacity = 0.7,
              group = "IMD",
              highlight = highlightOptions(
                dashArray = "",
                fillOpacity = 0.8,
                weight = 3,
                color = "Grey",
                bringToFront = TRUE)) %>%
  
  addPolygons(fillColor = ~pal4(LE.females.2015.2017),
              color = "white",
              weight = 2,
              opacity = 1,
              dashArray = "3",
              popup = popFLE,
              fillOpacity = 0.7,
              group = "FLE",
              highlight = highlightOptions(
                dashArray = "",
                fillOpacity = 0.8,
                weight = 3,
                color = "Grey",
                bringToFront = TRUE)) %>%
  
  addPolygons(fillColor = ~pal5(LE.males.2015.2017),
              color = "white",
              weight = 2,
              opacity = 1,
              dashArray = "3",
              popup = popMLE,
              fillOpacity = 0.7,
              group = "MLE",
              highlight = highlightOptions(
                dashArray = "",
                fillOpacity = 0.8,
                weight = 3,
                color = "Grey",
                bringToFront = TRUE)) %>%
  
  addPolygons(fillColor = ~pal6(GCSE.2018.A..C..),
              color = "white",
              weight = 2,
              opacity = 1,
              dashArray = "3",
              popup = popGCSE,
              fillOpacity = 0.7,
              group = "GCSE",
              highlight = highlightOptions(
                dashArray = "",
                fillOpacity = 0.8,
                weight = 3,
                color = "Grey",
                bringToFront = TRUE
              )) %>%
  
  addPolygons(fillColor = ~pal7(Leave),
              color = "white",
              weight = 2,
              opacity = 1,
              dashArray = "3",
              popup = popBrexit,
              fillOpacity = 0.7,
              group = "Brexit",
              highlight = highlightOptions(
                dashArray = "",
                fillOpacity = 0.8,
                weight = 3,
                color = "Grey",
                bringToFront = TRUE)) %>%
  
  #adding legends for each polygon
  addLegend(
    #setting the pallete to predetermined one
    pal=pal1,
    #where to get the values from        
    values = UKNUTS3MapWGS$GVA.in..2017,
    #group this with the GVA polygon        
    group = c("GVA"),
    #setting the title of the legend        
    title = "GVA in 2017 as a percentage of the UK average",
    #setting the position as bottom left      
    position = "bottomleft",
    #formatting the label so that it only shows to three significant figures to make it look cleaner     
    labFormat = labelFormat(digits=3))%>%
  
  #this is then replicated for the other indicators
  addLegend(pal=pal2,
            values = as.numeric(UKNUTS3MapWGS$U.E.rate...Jul.2018.Jun.2019),
            group = c("UE"),
            title = "Unemployment rate",
            position = "bottomleft",
            labFormat = labelFormat(digits=3))%>%
  
  addLegend(pal=pal3,
            values = UKNUTS3MapWGS$IMD.2019...Average.score,
            group = "IMD",
            title = "Indices of multiple deprivation score",
            position = "bottomleft",
            labFormat = labelFormat(digits=3))%>%
  
  addLegend(pal=pal4,
            values = UKNUTS3MapWGS$LE.females.2015.2017,
            group = "FLE",
            title = "Female life expectancy(years)",
            position = "bottomleft",
            labFormat = labelFormat(digits=3))%>%
  
  addLegend(pal=pal5,
            values = UKNUTS3MapWGS$LE.males.2015.2017,
            group = "MLE",
            title = "Male life expectancy(years)",
            position = "bottomleft",
            labFormat = labelFormat(digits=3))%>%
  
  addLegend(pal=pal6,
            values = UKNUTS3MapWGS$GCSE.2018.A..C..,
            group = "GCSE",
            title = "Percentage of 5 A*-C grades at GCSE",
            position = "bottomleft",
            labFormat = labelFormat(digits=3))%>%
  
  addLegend(pal=pal7,
           values = UKNUTS3MapWGS$Leave,
           group = "Brexit",
           title = "Percentage voting leave in 2016 EU referendum",
           position = "bottomleft",
           labFormat = labelFormat(digits=3))%>%
  
  #adding the option to control the layers
  addLayersControl(
    #the only base layer specified is the open streetmap layer
    baseGroups = "OSM (default)",
    #calling the overlay groups defined befined
    overlayGroup = c("UE", "GVA", "IMD", "FLE", "MLE", "GCSE", "Brexit"),
    #the collapsed = FALSE means that the options are always there
    options = layersControlOptions(collapsed = FALSE))%>%
  #initially hide every indicator but GVA
  hideGroup(c("UE", "IMD", "FLE", "MLE", "GCSE", "Brexit"))
              
#outputting the map
map

These plots are only a visual representation of the North-South divide compared to the three divides envisioned. Therefore, to better understand the divide we can use regression models. This is done using a dummy variable of 1 as being North and 0 being South according to the difference conceptions of the North-South divide. While many other factors are likely to influence these regression outputs, by using only a dummy variable we are able to explore to what extent can being in either the North or the South be associated with differences in indicators that can be broadly encompassed within measures related to quality of life. This is also able to comment on how much a simple dummy variable of North or South can capture the variance in these indicators through looking at the the R-squared value.

library(broom)

#setting up a function to run a linear regression
regressionmodel <- function(data1,data2){
  #setting the regression function
  modelx <- lm(data1 ~ data2)
  #assigning the summary to a variable
  modelx_summary <- summary(modelx)
  modelx_res <- tidy(modelx)
  #returning the summary of the regression
  return(modelx_summary)
}

#this regression model is then used on the different indicators according the divides of Dorling, Rowthorn and the South-East

#the first concern is the GVA outliers
hist(NUTS3_data$GVA.in..2017,
     main = "Histogram of GVA in 2017",
     xlab = "GVA",
     ylab = "frequency",
     col = "green",
     breaks = 100)

#this suggests there are outliers of four areas above 200% of GVA
#so these can be removed when running the regression

GVA_no_outliers <- subset.data.frame(NUTS3_data, NUTS3_data$GVA.in..2017<200)

GVA_outliers <- subset.data.frame(NUTS3_data, NUTS3_data$GVA.in..2017>200)

library(broom)

#GVA without extremes
GVA_Dorling <- regressionmodel(GVA_no_outliers$GVA.in..2017, GVA_no_outliers$Dorling)
GVA_Dorling
## 
## Call:
## lm(formula = data1 ~ data2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -36.349 -17.249  -5.249   9.751  77.951 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   95.949      2.645  36.281  < 2e-16 ***
## data2        -14.424      3.978  -3.626 0.000416 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 22.44 on 127 degrees of freedom
## Multiple R-squared:  0.09379,    Adjusted R-squared:  0.08666 
## F-statistic: 13.14 on 1 and 127 DF,  p-value: 0.000416
GVA_Rowthorn <- regressionmodel(GVA_no_outliers$GVA.in..2017, GVA_no_outliers$Rowthorn)
GVA_Rowthorn
## 
## Call:
## lm(formula = data1 ~ data2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -34.975 -18.171  -4.071   9.825  79.325 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   94.575      2.594  36.461  < 2e-16 ***
## data2        -12.404      4.085  -3.036  0.00291 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 22.76 on 127 degrees of freedom
## Multiple R-squared:  0.06767,    Adjusted R-squared:  0.06033 
## F-statistic: 9.218 on 1 and 127 DF,  p-value: 0.002908
GVA_SE <- regressionmodel(GVA_no_outliers$GVA.in..2017, GVA_no_outliers$SE)
GVA_SE
## 
## Call:
## lm(formula = data1 ~ data2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -35.998 -12.842  -3.698   9.402  73.202 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  100.698      2.864  35.162  < 2e-16 ***
## data2        -19.656      3.807  -5.163 9.12e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 21.43 on 127 degrees of freedom
## Multiple R-squared:  0.1735, Adjusted R-squared:  0.167 
## F-statistic: 26.66 on 1 and 127 DF,  p-value: 9.121e-07
#GVA with extremes
GVA_Dorling_E <- regressionmodel(NUTS3_data$GVA.in..2017, NUTS3_data$Dorling)
GVA_Dorling_E
## 
## Call:
## lm(formula = data1 ~ data2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
##  -67.93  -38.53  -16.62    0.88 1109.67 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   127.53      14.41    8.85 5.23e-15 ***
## data2         -46.00      22.01   -2.09   0.0386 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 125.6 on 131 degrees of freedom
## Multiple R-squared:  0.03227,    Adjusted R-squared:  0.02488 
## F-statistic: 4.368 on 1 and 131 DF,  p-value: 0.03856
GVA_Rowthorn_E <- regressionmodel(NUTS3_data$GVA.in..2017, NUTS3_data$Rowthorn)
GVA_Rowthorn_E
## 
## Call:
## lm(formula = data1 ~ data2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
##  -64.67  -39.97  -18.17    0.93 1112.93 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   124.27      14.00   8.876 4.52e-15 ***
## data2         -42.10      22.39  -1.880   0.0623 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 126 on 131 degrees of freedom
## Multiple R-squared:  0.02628,    Adjusted R-squared:  0.01884 
## F-statistic: 3.535 on 1 and 131 DF,  p-value: 0.0623
GVA_SE_E <- regressionmodel(NUTS3_data$GVA.in..2017, NUTS3_data$SE)
GVA_SE_E
## 
## Call:
## lm(formula = data1 ~ data2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
##  -75.68  -35.98  -11.44    2.76 1096.82 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   140.38      16.03   8.756 8.84e-15 ***
## data2         -59.34      21.64  -2.742  0.00696 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 124.2 on 131 degrees of freedom
## Multiple R-squared:  0.05428,    Adjusted R-squared:  0.04706 
## F-statistic: 7.519 on 1 and 131 DF,  p-value: 0.00696
#unemployment
UE_Dorling <- regressionmodel(NUTS3_data$U.E.rate...Jul.2018.Jun.2019, NUTS3_data$Dorling)
UE_Dorling
## 
## Call:
## lm(formula = data1 ~ data2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.7255 -0.8497 -0.0056  0.7253  3.5745 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   3.8222     0.1352  28.277  < 2e-16 ***
## data2         0.7033     0.2044   3.441 0.000786 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.147 on 126 degrees of freedom
##   (5 observations deleted due to missingness)
## Multiple R-squared:  0.08591,    Adjusted R-squared:  0.07866 
## F-statistic: 11.84 on 1 and 126 DF,  p-value: 0.0007861
UE_Rowthorn <- regressionmodel(NUTS3_data$U.E.rate...Jul.2018.Jun.2019, NUTS3_data$Rowthorn)
UE_Rowthorn
## 
## Call:
## lm(formula = data1 ~ data2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.6566 -0.8566 -0.0517  0.6865  3.6434 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   3.9135     0.1333  29.367   <2e-16 ***
## data2         0.5431     0.2111   2.572   0.0113 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.169 on 126 degrees of freedom
##   (5 observations deleted due to missingness)
## Multiple R-squared:  0.04989,    Adjusted R-squared:  0.04235 
## F-statistic: 6.617 on 1 and 126 DF,  p-value: 0.01126
UE_SE <- regressionmodel(NUTS3_data$U.E.rate...Jul.2018.Jun.2019, NUTS3_data$SE)
UE_SE
## 
## Call:
## lm(formula = data1 ~ data2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.5149 -0.8188 -0.0149  0.6851  3.7851 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   3.9202     0.1527  25.666   <2e-16 ***
## data2         0.3947     0.2096   1.883    0.062 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.183 on 126 degrees of freedom
##   (5 observations deleted due to missingness)
## Multiple R-squared:  0.02738,    Adjusted R-squared:  0.01966 
## F-statistic: 3.547 on 1 and 126 DF,  p-value: 0.06196
#Male life expectancy
MaleLE_Dorling <- regressionmodel(NUTS3_data$LE.males.2015.2017, NUTS3_data$Dorling)
MaleLE_Dorling
## 
## Call:
## lm(formula = data1 ~ data2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.2509 -0.8518  0.0473  0.7982  2.4491 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  80.0527     0.1366 586.230  < 2e-16 ***
## data2        -1.6018     0.2070  -7.738 2.56e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.175 on 129 degrees of freedom
##   (2 observations deleted due to missingness)
## Multiple R-squared:  0.317,  Adjusted R-squared:  0.3117 
## F-statistic: 59.87 on 1 and 129 DF,  p-value: 2.563e-12
MaleLE_Rowthorn <- regressionmodel(NUTS3_data$LE.males.2015.2017, NUTS3_data$Rowthorn)
MaleLE_Rowthorn
## 
## Call:
## lm(formula = data1 ~ data2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.2154 -0.8450  0.0253  0.7550  2.4846 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  79.9747     0.1345 594.576  < 2e-16 ***
## data2        -1.5593     0.2135  -7.304 2.57e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.196 on 129 degrees of freedom
##   (2 observations deleted due to missingness)
## Multiple R-squared:  0.2926, Adjusted R-squared:  0.2871 
## F-statistic: 53.35 on 1 and 129 DF,  p-value: 2.572e-11
MaleLE_SE <- regressionmodel(NUTS3_data$LE.males.2015.2017, NUTS3_data$SE)
MaleLE_SE
## 
## Call:
## lm(formula = data1 ~ data2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.5740 -0.9379  0.0260  1.1190  2.2260 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  80.0879     0.1655 484.047  < 2e-16 ***
## data2        -1.3140     0.2216  -5.928 2.63e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.26 on 129 degrees of freedom
##   (2 observations deleted due to missingness)
## Multiple R-squared:  0.2141, Adjusted R-squared:  0.208 
## F-statistic: 35.14 on 1 and 129 DF,  p-value: 2.63e-08
#Female life expectancy
femaleLE_Dorling <- regressionmodel(NUTS3_data$LE.females.2015.2017, NUTS3_data$Dorling)
femaleLE_Dorling
## 
## Call:
## lm(formula = data1 ~ data2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.55549 -0.74644  0.05903  0.63056  2.78976 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  83.7562     0.1201 697.328  < 2e-16 ***
## data2        -1.6809     0.1835  -9.161 9.06e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.047 on 131 degrees of freedom
## Multiple R-squared:  0.3905, Adjusted R-squared:  0.3859 
## F-statistic: 83.93 on 1 and 131 DF,  p-value: 9.061e-16
femaleLE_Rowthorn <- regressionmodel(NUTS3_data$LE.females.2015.2017, NUTS3_data$Rowthorn)
femaleLE_Rowthorn
## 
## Call:
## lm(formula = data1 ~ data2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.58406 -0.73617  0.05051  0.64223  2.88603 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  83.6599     0.1208 692.285  < 2e-16 ***
## data2        -1.5963     0.1933  -8.259 1.39e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.088 on 131 degrees of freedom
## Multiple R-squared:  0.3424, Adjusted R-squared:  0.3374 
## F-statistic: 68.22 on 1 and 131 DF,  p-value: 1.395e-13
femaleLE_SE <- regressionmodel(NUTS3_data$LE.females.2015.2017, NUTS3_data$SE)
femaleLE_SE
## 
## Call:
## lm(formula = data1 ~ data2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.87214 -0.80257 -0.00812  0.83224  2.72678 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  83.8192     0.1464  572.35  < 2e-16 ***
## data2        -1.4272     0.1977   -7.22 3.79e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.134 on 131 degrees of freedom
## Multiple R-squared:  0.2847, Adjusted R-squared:  0.2792 
## F-statistic: 52.13 on 1 and 131 DF,  p-value: 3.786e-11
#GCSE
GCSE_Dorling <- regressionmodel(NUTS3_data$GCSE.2018.A..C.., NUTS3_data$Dorling)
GCSE_Dorling
## 
## Call:
## lm(formula = data1 ~ data2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15.9579  -2.9579   0.1421   3.5421  11.7039 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  64.9961     0.5822 111.642  < 2e-16 ***
## data2        -3.4382     0.8893  -3.866 0.000173 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.075 on 131 degrees of freedom
## Multiple R-squared:  0.1024, Adjusted R-squared:  0.09556 
## F-statistic: 14.95 on 1 and 131 DF,  p-value: 0.0001733
GCSE_Rowthorn <- regressionmodel(NUTS3_data$GCSE.2018.A..C.., NUTS3_data$Rowthorn)
GCSE_Rowthorn
## 
## Call:
## lm(formula = data1 ~ data2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -19.0654  -2.8654  -0.0423   3.2346  12.0346 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  64.6654     0.5734 112.774  < 2e-16 ***
## data2        -2.9231     0.9170  -3.188  0.00179 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.161 on 131 degrees of freedom
## Multiple R-squared:  0.07198,    Adjusted R-squared:  0.06489 
## F-statistic: 10.16 on 1 and 131 DF,  p-value: 0.001794
GCSE_SE <- regressionmodel(NUTS3_data$GCSE.2018.A..C.., NUTS3_data$SE)
GCSE_SE
## 
## Call:
## lm(formula = data1 ~ data2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -16.4274  -3.0417   0.5726   3.8726  11.3583 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  65.3417     0.6575  99.382  < 2e-16 ***
## data2        -3.3143     0.8875  -3.735 0.000279 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.093 on 131 degrees of freedom
## Multiple R-squared:  0.09622,    Adjusted R-squared:  0.08932 
## F-statistic: 13.95 on 1 and 131 DF,  p-value: 0.0002795
#IMD 
IMD_Dorling <- regressionmodel(NUTS3_data$IMD.2019...Average.score, NUTS3_data$Dorling)
IMD_Dorling
## 
## Call:
## lm(formula = data1 ~ data2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -14.6288  -4.4745  -0.1967   4.2287  18.6832 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  19.0696     0.7748  24.613  < 2e-16 ***
## data2         7.2862     1.1805   6.172 8.22e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.665 on 128 degrees of freedom
##   (3 observations deleted due to missingness)
## Multiple R-squared:  0.2294, Adjusted R-squared:  0.2233 
## F-statistic:  38.1 on 1 and 128 DF,  p-value: 8.223e-09
IMD_Rowthorn <- regressionmodel(NUTS3_data$IMD.2019...Average.score, NUTS3_data$Rowthorn)
IMD_Rowthorn
## 
## Call:
## lm(formula = data1 ~ data2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -14.8630  -4.7845  -0.0544   4.4738  18.4490 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  19.3796     0.7552   25.66  < 2e-16 ***
## data2         7.2104     1.2057    5.98 2.08e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.712 on 128 degrees of freedom
##   (3 observations deleted due to missingness)
## Multiple R-squared:  0.2184, Adjusted R-squared:  0.2123 
## F-statistic: 35.76 on 1 and 128 DF,  p-value: 2.081e-08
IMD_SE <- regressionmodel(NUTS3_data$IMD.2019...Average.score, NUTS3_data$SE)
IMD_SE
## 
## Call:
## lm(formula = data1 ~ data2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -13.1443  -6.1962   0.2008   4.9998  20.1713 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  19.1056     0.9061  21.085  < 2e-16 ***
## data2         5.7621     1.2348   4.666 7.62e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.019 on 128 degrees of freedom
##   (3 observations deleted due to missingness)
## Multiple R-squared:  0.1454, Adjusted R-squared:  0.1387 
## F-statistic: 21.77 on 1 and 128 DF,  p-value: 7.616e-06
#Brexit leave
Brexit_Dorling <- regressionmodel(NUTS3_data$Leave, NUTS3_data$Dorling)
Brexit_Dorling
## 
## Call:
## lm(formula = data1 ~ data2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -29.225  -4.393   0.575   6.707  21.675 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   50.625      1.112  45.510  < 2e-16 ***
## data2          6.868      1.699   4.042 8.99e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.698 on 131 degrees of freedom
## Multiple R-squared:  0.1109, Adjusted R-squared:  0.1041 
## F-statistic: 16.34 on 1 and 131 DF,  p-value: 8.992e-05
Brexit_Rowthorn <- regressionmodel(NUTS3_data$Leave, NUTS3_data$Rowthorn)
Brexit_Rowthorn
## 
## Call:
## lm(formula = data1 ~ data2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -29.7086  -4.1000   0.2914   6.5914  21.1914 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   51.109      1.090   46.90  < 2e-16 ***
## data2          6.291      1.743    3.61 0.000435 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.808 on 131 degrees of freedom
## Multiple R-squared:  0.09047,    Adjusted R-squared:  0.08352 
## F-statistic: 13.03 on 1 and 131 DF,  p-value: 0.0004352
Brexit_SE <- regressionmodel(NUTS3_data$Leave, NUTS3_data$SE)
Brexit_SE
## 
## Call:
## lm(formula = data1 ~ data2)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -27.98  -5.68   0.32   6.82  22.92 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   49.380      1.232  40.067  < 2e-16 ***
## data2          7.631      1.664   4.587 1.04e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.546 on 131 degrees of freedom
## Multiple R-squared:  0.1384, Adjusted R-squared:  0.1318 
## F-statistic: 21.04 on 1 and 131 DF,  p-value: 1.038e-05

The results in these regressions were then inputted into figure 4 in the PDf.